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' Abstract 
Or 

A probability distribution function is presented which provides a realistic description of the detection of 

m : 

scattered photons. The resulting probabilities can be described analytically by means of a superposition 
' of several special functions. These exact expressions can be evaluated numerically only for small 

distances and limited time residuals, due to computer accuracy limitations. In this report we provide 
approximations for the exact expressions in different regions of the distance-time residual space, defined 
by the detector geometry and the space-time scale of an event. These approximations can be evaluated 
numerically with a relative error with respect to the exact expression at the boundaries of less than 

10- 3 . 
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1 Introduction 



■ In the track reconstruction process in neutrino telescopes it is customary to perform a 

log-likelihood minimization using a gamma function based Probability Distribution Function 
(PDF). The use of a so-called Pandel PDF [1] has been investigated in the reconstruction 
of muon tracks in the AMANDA detector [2]. However, the problem with a Pandel PDF is 
^ . that it cannot cope with negative time residuals which may result due to time jitter in the 

detection devices. One way to overcome this limitation is to consider a convolution of a Pandel 
PDF with a Gaussian, the latter accounting for the finite time resolution of the detector. The 
resulting expression, referred throughout this report as a CPandel PDF, can be evaluated 
analytically. This involves the superposition of various special functions, which are available 
in the GNU Scientific Library (GSL) [3]. However, due to computer accuracy limitations, the 
numerical evaluation of a CPandel PDF is restricted to a rather limited distance-time residual 
domain. In this report we present analytical approximations which are necessary to extend 
the applicability domain of a CPandel PDF for practical cases. The implementation shown 
here is tailored for the geometry of the AMANDA neutrino telescope [2] and has been realized 
in an analysis framework (IcePack) based on ROOT [4] and an analysis toolbox (Ralice) [5] 
which was originally developed for the Alice experiment at the future CERN-LHC collider. 
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2 The CPandel probability distribution function 



In a reconstruction procedure, track parameters emerge as the solution of an optimization 
problem: given experimentally measured values (in our case photon detection times thu), find 
the values of parameters minimizing the log of the likelihood C. The latter is given by [2] 

£ = II P( a > thitj) , (!) 
3 

where a denotes the parameter(s) characterizing the hypothesis of an event and the PDF 
p(ci, thitj) describes a photon arrival at the receiver j. It is convenient to use the time residual 
t as a PDF variable 

t = ifiit tgeom i (2) 

where t geom is the photon arrival time for a case of no scattering and no absorption. 

Note that in our notation we have dropped the index j, since for the scope of the present 

report it is sufficient to limit ourselves to a one source-one receiver situation. 

The generation of detector signals corresponds to a process in which a charged particle 
moves along a straight track with velocity v exceeding the in-medium light speed c m . This 
gives rise to emission of Cherenkov radiation, which is recorded by the detection devices. The 
geometrical situation is sketched in Fig. 1, in which the particle starts from a point X at 
(r*o,to) an d arrives at a point B when the Cherenkov front hits the detector at (rhit,thit)- 
The track point labeled E indicates the point of closest approach w.r.t. the receiver location. 
The geometrical (expected) arrival time can be evaluated by considering Cherenkov photons 
emitted from the particle track and propagating freely (without scattering and absorption) 
to arrive at the receiver at (r hit, tut)- This specific case implies that t geom = thit- 




EB :At v EG : At -v . EF : At -v 

^ phase group 

Fig. 1. Geometry of the signal generation process. Further details can be found in the text. 

Indicating by At the time it takes for the particle to travel from E to B, it is seen from 
Fig. 1 that in the same time interval the Cherenkov front has moved with phase velocity 
Vphase = c/n p h from E to G. Here c indicates the light speed in vacuum and n p h is called the 
phase refractive index. However, a detector is sensitive for real photons which travel with the 



2 



group velocity v 



group 



c/rig r , where n gr is called the group refractive index. 



Since v group < v p hase, the real photon signal lags behind the Cherenkov front. In Fig. 1 this is 
indicated by the line element EF, traveled by the real Cherenkov photons in the time interval 
At mentioned above. So, one can regard the real photons to be located on a wavefront given by 
BF which is comparable to the Cherenkov front provided the complement of the Cherenkov 
angle 9 C is reduced by a. 

Assuming v = c we have cos(6> c ) = l/n p h. The expected arrival time is then given by 



geom 



t + 



1 



-, . , WgrTlph 1 

v ■ r + a , — 



(3) 



where v indicates the unit vector in the moving direction of the particle. 



Minimization w.r.t. the time residuals t = 



L geom 



for all observed signals provides the 



basis for track reconstruction, as mentioned before. 

To achieve this, a Pandel PDF p, t) was suggested in [1] : 



no 



(4) 



where £ and p are phenomenological parameters related to the characteristics of the medium. 
The parameter £ represents the distance between the emission and detection locations of 
a Cherenkov photon in units of the mean photon scattering length A, i.e. £ = d/(Asin# c ). 
Based on the experiences from the muon track reconstruction procedure with the AMANDA 
Neutrino Telescope [2], we use the following parameter values throughout this report 



A = 33.3 m, p = 0.004 ns' 



(5) 



A PDF for realistic signals should account for the finite time resolution of the detector. 
This can be achieved by convolving a Pandel PDF p (which describes only the propagation 
of photons in a medium, i.e. an ideal measurement situation of a detector with a zero time 
resolution) with a time jitter function. The latter may be described by a gaussian with a mean 
of zero and standard deviation a. The width of the gaussian distribution, a, represents the 
time resolution of the detector and usually is made up from various sources, which motivates 
the use of a gaussian. We call this PDF CPandel and denote it as T a : 



dx 

V / 27T(T 2 



(6) 



The convoluted PDF T a (p,^,t) can be calculated exactly [6] : 



2( 1 +0/2 



r(*fc + i)) 



(7) 



where 



rj = pa 



(7 



(8) 



and \F\ is the confluent hypergeometric function [7]. The latter is implemented in the GNU 
Scientific Library (GSL) [3], which enables numerical evaluation of J : a (p^,t). 
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Accounting for the detector finite time resolution provides a solution of the problem of 
negative time residuals : T a exists for any t and causality is satisfied. For t < the finite 
value of the CPandel PDF is purely due to the detector time jitter [6]. A typical profile of T a 
is shown in Fig. 2. 




-30 -20 20 40 60 80 100 120 

t(ns) 



Fig. 2. CPandel PDF for two values of the time resolution parameter a (5 and 10 ns). The distance of 
closest approach d is 10 m. 

In a reconstruction procedure the minimizer moves randomly within the £ — t (i.e. distance- 
time residual) plane. As such, it is essential to maximize the area in which numerical evalua- 
tions of T a can be performed. 

Although T a possesses all the desirable attributes for a PDF [6] and \F\ is implemented in 
GSL, computer accuracy limitations pose a problem in minimizing the log- likelihood based on 
T a - For instance, the attempt to calculate .7-5(0. 004, 2, 190) fails due to overflow in the GSL 
procedure. Consequently, the GSL implementation of \F\ does not allow evaluation of T a 
when the distance is 66.6 m and the time residual is 190 ns. As such, evaluation of the exact 
expression of T a as given by (7) cannot be performed in all the areas relevant for realistic 
cases. 

Due to these limitations, it is necessary to find an approximation which allows numerical 
evaluation of the CPandel PDF in all physically relevant areas, while retaining all the qual- 
ities defined from the exact expression (7). Our strategy is to use (7) in the area where the 
GSL implementation of \F\ yields reliable results and outside that area use an approximation 
which allows extension of the support for the minimizer procedure. 

Since the analytic properties of \F\ are well established [7], it is straightforward to find 
approximations for T a in terms of uniform (i.e. numerically convergent) expansions 0- 
Depending on the position in the £ — t plane (as the minimizer moves away from the origin), 
different functions have to be used to approximate T a - However, it should be noted that these 
different expressions for the CPandel PDF (corresponding to different regions in the distance- 
time plane) represent the same function approximated by means of analytic continuation and 
uniform expansion. In other words, the results presented hereafter are not originating from 
introducing any new ad hoc functions or parameters. 

1 Note that using an asymptotic expansion (with regard to a variable, say t) for the CPandel PDF 
may disturb the approximation when another variable (£) is varied [6] . 
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Below we present the various expressions for the CPandel PDF corresponding to different 
relevant regions of the distance-time residual plane. 



3 Expressions for T a in different distance-time regions 

For direct hits (i.e. small £ and small |t|), we can use the exact expression for JF CT as given 
by eq. (7), since this is correctly handled by the GSL implementation on this domain. 
To extend the coverage to all relevant areas in the £ — t plane, it is necessary and sufficient to 
consider the following regimes for the variables £ and n (for n see (8)) : positive and negative r\ 
(corresponding to t < pa 2 and t > pa 2 , respectively), £ < 1 and £ > 1. Approximations for T a 
in the latter four regions exhaust all the possibilities following from the analytic continuation 
procedure of the function \F\, of which the details can be found at the end of this report. 
Consequently, the support for T a consists of the following five parts : 

Region 1 : direct hits region - small £ and small \t\. 
Region 2 : large positive t and restricted £. 
Region 3 : large £ and rj < (i.e. t > per 2 ). 
Region 4 : large £ and r] > (i.e. t < pa 2 ). 
Region 5 : large negative t and restricted £. 

In the following we denote the CPandel PDF in region j by fj. 

Zero distance 

On the i-axis, i.e. where £ = 0, the exact expression (7) turns into 



which we use for £ = 0. 

Region 1 : small £, small \t\ 

Here GSL allows the use of the exact expression (7). In other words, in the region corre- 
sponding to direct hits no approximation is used. 
In this region, for t < we use the constraint 



which yields satisfactory results in matching the values of T a calculated via the exact expres- 
sion (7) with the values of f± and fa (as outlined below). 

Region 2 : £ < 1 and large positive t (i.e. t > a) 

In this region, T a is approximated by 




(9) 



V2ira 2 



-5a<t<0 



(10) 



Fa(p,Z,t) = e^ 2 / 2 p(p,tt) = / 2 (p,£,t) • 



(11) 



Here p(p, £,i) is the Pandel PDF as given in eq. (4). 
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Region 3 : £ > 1 and t > pa 2 

In this region T a is approximated by 



where 



(12) 



v 2 £ i 



1 



<> --^2 + t - I + 4 + fc(2e - 1} " 4 ln(1 + " 2) " I ln(2) 

£-1 



+ ^- ln(2£ - 1) + £ln(p) + (£ - 1) ln(<7) , 



Vl + z 2 + ln^z + v 7 ! + 



> o , 



$ = 1 - 



m n 2 



(2£-l) (2£-l) 2 ' 
288 



JV X = A (20/3 2 + 30/3 + 9), iV 2 = ^ (6160/3 4 + 18480/3 3 + 19404/3 2 + 8028 / + 9 15) 



1 



P = o 



2 vvtt^ 2 



- 1 . 



Region 4 : £ > 1 and t < pa 2 

In this region, T a is approximated by 



2vr 



[/(£) e~ fc C*" 1 ) (1 + z 2 )" 1 / 4 * = / 4 (p, £, t) , (13) 



where 



[/(£) = e^ 2 - 1 / 4 (2£ - l)-^ 2 2«- 4 )/ 2 , 



z Vl + z 2 + In + v 7 ! + z 2 



iVi iV 2 
* = 1 + 7X7 rr + 



(2£-l) (2£-l) 2 • 
Nj , and f3 are the same as for region 3 and z is now defined as 

> . 



V4£^2 



Region 5 : £ < 1 and t <C —a 



In this region T a is approximated by 



Fa(p,£,t) 



V2-KO- 2 



(14) 
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4 Description of the support of T a in the £ — t plane 



Concrete limiting values of the distances and time residuals denning the above five regions 
depend on the values of the Pandel PDF parameters and on the value of the time jitter 
parameter a. Values of A and p are given by (5) and according to realistic time jitter values 
[2], we consider the cases a = 5 and a = 10 ns. 

The corresponding areas of the £ — t plane where T a can be evaluated are presented in Figs. 3 
and 4. 
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Fig. 3. Support of T a for a = 5 ns. Distances and time residuals are restricted to the rectangle 
< £ < 30 (0 < d < 999 m) and -150 < t < 3500 ns (pa 2 = 0.1 ns). 



Outside the large rectangle, numerical evaluation of T a becomes impossible. With increas- 
ing distance, T a decreases sharply to values of about 10~ 200 — 10 -280 , depending on the value 
of t. For fixed £ and increasing \t\, computer limitations lead to overflow. 
To compensate for these effects, one can evaluate T a on the corresponding boundary location 
and impose a penalty value on the resulting logdikelihood value beyond the regions 1 — 5. 

The accuracy while crossing the border of the regions i and j is given by inequality: 

It + \ — t) ~ t)\ s 1n -3 /ir\ 

r ij,k{£, t) = < 1U , (15) 

where k = i or k = j. 

The value of the relative error of about 10~ 3 occurs while comparing the function f\ with fa 
and /4 with /s, i.e. the region with negative t and small £. For the remaining areas, the value 
of rjj is in the range of 10~ n < r.- L j < 10 -6 . 
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Fig. 4. Support of T„ for a — 10 ns. Distances and time residuals are restricted to the rectangle 
0<£<50(0<d< 1665 m) and -250 < t < 3500 ns {pa 2 = 0.4 ns). 

5 Derivations of the various approximations for J- a (p,^,t) 

Approximation in the region 5: £ < 1 and t <C —a 

In this region r/ = pa — t/a is positive and we can use the asymptotic expansion formula 
13.5.1 (bold faced numeration refers to the corresponding sections of [7]) valid for large z 
and fixed a, b : 



iF!(a,b,z) e l7Ta z- a 



T(b) 



T(b-a) 



R-l 



T{a + n)T(l + a-b + n) 



^ T(a)T(l + a-b)n\ 



~ n + 0{z 



+ 



e z z a-b 



5-1 

E 

n=0 



T(b-a + n)T(l-a + n) 



T(b-a)T(l 



ami 



+ 0{z 



(16) 



Application of (16) to the combination of hypergeometric functions appearing in the exact 
expression (7) for Jv(p, £,i) : 

r(|(e + i)) v v nho { ] 

shows that the rising terms, proportional to e z , originating from the second term of (16), 
cancel. For the leading order term we obtain 



r(i)^- 1 2( 1 -g)/ 2 r(i)^- 1 2( 1 -^/ 2 
rxJorxW+i)) ~r(k)ra(e + i)) 



o 



(18) 
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and the same cancellation occurs for the non-leading order terms. 

The remaining terms of (16), proportional to e lna , give rise (in leading order) to the ex- 
pression 



r(i(e + i))r(i(i-0) t{\oy{\{z + i)) 

1 e™' 2 

r(i(e + i))r(i(i-0) T{\onW + V) 



(19) 



Using the reflection formulae 6.1.17 : 



r(z)r(i - z )= j . and r(z + -)r(z - -) = — ^— , (20) 

v ' v ; sin(Trz) v 2' v 2 ; cos(vrz) ' v 7 

we can write the leading order expression (19) as 

-!■ I cos«/2) - e™l 2 sin«/2) = ^ . (21) 



7T \ / 7T 

So, the imaginary part is vanishing, as is expected for a PDF. 

Using (21) in the exact expression (7) for F a (p., £, t) results in the leading order approximation 
(14) used in the region 5. 

Approximation in the region 2: £ < 1 and t ^> a 

Also here we derive only the leading order approximation using (16). 
When t > pa 2 , n is negative and for convenience we introduce a positive p : 

fi = -r) = --pa>0 . (22) 
a 

The approximation in the region 2 is obtained by applying (16) to the exact expression (7) 
for J~a(Pi£,,t), written in terms of this positive p. The net effect is that now the sign in (17) 
between the two i.Fi's is reflected. Lengthy but straightforward algebra leads to 

„« ^-i e p^n-,i 2 i-s r(i) ft - p* 2 ^- 1 



Using the duplication formula 6.1.18 : 



T(2z) = 2 2z - 1 tt- 1 / 2 T(z)T(z + i) (24) 



we arrive at the approximation in the region 2 : 



^(M,t)^ 2 ^^e^ 



(25) 



In deriving the approximations for the regions 2 and 5 the non-leading terms, proportional 
to (pa 2 /t)~ n , have been neglected. The corresponding numerical values contributing to the 
correction factors amount to less than 0.1%. Note that the asymptotic expansion formula (16) 
allows for the accounting of the non-leading terms in a straightforward way if the necessity 
arises. 
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Approximations in the regions 3 and 4 

To proceed further it is convenient to express the PDF J- a {p, £, t) in terms of the parabolic 
cylinder function U(a, x) [7]. Comparison with 13.6.15 and 13.6.16 shows that the PDF (7) 
can be expressed as 



(pa)* 

exp 



V2 



7T(T Z 



e_ 

2a 2 ) 



(26) 



However, the function U is not implemented in GSL, therefore we can not use the above 
expression in the region of small £ and small ±t. This disadvantage is compensated by the 
fact that manipulation with U leads to more compact expressions than those with the su- 
perposition of the two hypergeometric functions. E.g., it is straightforward to verify that the 
asymptotic expansion formula for U, 19.8.1, immediately leads to the results for the regions 
2 and 5, and there is no need to follow up the intermediate steps of calculation in order to 
see explicitly that the rising terms cancel as we saw in (18), or that the result contains no 
imaginary part, as we saw in (21). 
Let us consider the region 4, where t < pa 2 , i.e. where 



7] = pa > . 

a 

Below we list the essential steps necessary to obtain the approximation. 
Introduction of a new variable z, 

= 71 

transforms the equation for a parabolic cylinder function 19.1.2 into 
d 2 



d22 ~(e-i,) = P«-i) 2 (i+, 2 )t/(«-i 

Next, we introduce a function U\ which is defined as 

U^(l + z^u(t-\,z) 
and define the variables </> and k as 



(f) = sinh 1 (z) , k = - sinh(</>) cosh(</>) + ln(^sinh(0) + cosh(0) 



Finally, we introduce a function ^ and variable (3 as 



* = ekre-D Ult p= i(tanh(0 - 1) = \ y^-^ 



Making use of the relations 
dk 



dk 



— = cosh(^) = \/ z 2 + 1 and — = 2 cosh 4 (0) 
d(p dp 



1 



(3 2 (1 + P) 2 

it is straightforward to verify that the function ^ satisfies the following equation 



d_ 

dp 



(2$ - 1) d* 20/3 2 + 20/? + 3 = 
4 df3 16 



(27) 



(28) 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 
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Expanding ^ in inverse powers of 2£ — 1, we obtain 



W " 1 + 2£ - 1 + (2£ - l) 2 + {6b) 

The first term of the expansion (35) amounts to one, because we have introduced the trans- 
formation u(£ - §, z) -»■ Ui -»■ 

The expansion coefficients aj(/?) of (35) are obtained by using the expansion (35) of * in 
the equation (34). Neglecting the terms 0((2£ — l) -1 ), i.e. retaining only constant terms, we 
obtain the following relation for ai(/3) : 

dai(P) 20/? 2 + 20/3 + 3 

"dF = 4 • (36) 

The corresponding solution is given by 

ai(/3) = ^(20/? 2 + 30/? + 9) , (37) 

which is exactly the numerator Ni in the expression for the approximation of the PDF in 
region 4. Analyzing in the same fashion the (2£ — 1) _1 terms, we obtain the expression for 
a,2(f3), which corresponds to N2 in the approximation expression of the PDF in region 4. 
The approximation (13) of the PDF for the region 4 is obtained by collecting the various 
terms together, transforming back to the parabolic cylinder function * — ► Ui — ► U(£ — \, z) 
and using (28) and (26). 

The approximation (12) in the region 3 is obtained using the same strategy as the one used 
for the region 4. The difference we need to account for is that now 77 = pa — t/a is negative 
which results in a slightly different series for \£. 



6 Summary 

Based on approximations given by eqs. (11)-(14) we have evaluated the values of a Gauss 
convoluted Pandel (CPandel) PDF in an area of the distance-time residual plane which covers 
basically all physically relevant parameters to perform track reconstruction in a neutrino tele- 
scope. The approximations are obtained by considering analytical continuations of the exact 
expression (7) for the CPandel PDF in different areas of the £ — t plane and, wherever nec- 
essary, expanding the result in a numerically convergent series. As such, our approximations 
are obtained without involving any new ad hoc functions or parameters. 

The concrete values of £ and t which define the borders of the regions where each of the 
approximations is applicable, are defined by values of the parameters of the original Pandel 
PDF. In this report we have used the values given in (5). Since the expressions (11)-(14) are 
derived analytically, it is a matter of straightforward interpolation to define support for a 
CPandel PDF for different parameter values. 
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